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Summary 

The integral form of the complete, unsteady, compressible, three-dimensional Navier-Stokes 
equations in the conservation form, cast in generalized coordinate system, are solved, numerically, 
to simulate the vortex breakdown phenomenon. The inviscid fluxes are spatially discretized using 
Roe’s upwind-biased flux -difference splitting scheme and the viscous fluxes are discretized using 
central differencing. Time integration is performed using a backward Euler ADI scheme. A full 
approximation multigrid is used to accelerate the convergence to steady state. 

1. Introduction 

The phenomenon of 'vortex bursting" or "vortex breakdown" was first reported by Peckham 
and Atkinson^ in 1957. They observed that under certain conditions, the vortex core shed by the 
leading edge of the Gothic wing they were testing would swell, eventually stagnating the flow along 
the vortex axis forming a bubble of recirculating flow. In the following years, similar observations 
were made by Elle* 21 ^, Werle^, and Lambourne and Bryer 151 . The classic photograph in reference 
5, reproduced here as Fig. 1-1, shows the sudden breakdown of the leading-edge vortices, once 
a critical angle of attack 3S reached, to either a symmetric bubble-type structure or an asymmetric 
spiral-type structure. Downstream of this structural change, the flow is usually highly turbulent 
and diffusive. 



Figure 1-1. Leading-edge Vortex Breakdown (Lambourne and Bryer 1961) 
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The effect of vortex breakdown on the aerodynamics of a wing is very important. Hummel 
and Srinivasan* 6 ' showed that, with increasing angle of attack, the lift, the drag and the pitching 
moment of a delta wing deteriorates abruptly as the vortex breakdown location moved upstream 
over the trailing edge of the wing. Peake and Tobak* 7 ' pointed out that vortex breakdown causes 
buffeting, unsteadiness and poor control. Modem fighter aircraft, like the F-18, that depend on 
vortex-induced lift for high-alpha maneuvers, experience significant loss of performance -loss of 
lift increase in drag and an increase in nose-up pitching moment- at certain attitudes, due to vortex 
breakdown. The fatigue life of some of the aircraft components, like the vertical stabilizers, are 
drastically reduced due to vibrations caused by vortex breakdown. Vortex breakdown is also 
known to cause undesirable effects in turbomachinery. 

Some uses of vortex breakdown have also been suggested. Vortex breakdown may be used to 
dissipate the strong trailing vortices generated by large aircraft. The zone of recirculation flow of 
a vortex that has undergone breakdown can be used as a flame holder in gas turbine combustion 
chambers. Fuel and air can undergo intense turbulent mixing in the recirculating zone and result in 
a stable and compact flame. This will help increase combustion efficiency, reduce the combustor 
size and control the pollutants* 8 '. However the importance of vortex breakdown to aircraft designed 
to do high-alpha maneuvers is the main motivation of NASA to understand this phenomenon. 

The phenomenon of vortex breakdown has been studied, both experimentally and theoretically, 
for over 30 years. The most important contribution to the experimental study of vortex breakdown 
was made by Harvey* 9 ', who isolated the vortex from the wing. He successfully simulated vortex 
breakdown in tube flow with swirl introduced upstream by a system of guide vanes. This set 
the stage for most of the "clean" experimental work that followed. Theoretically, there had been 
many efforts to explain the bubble-type vortex breakdown. Although many theories have been 
proposed, none of them explain, satisfactorily, all the details of breakdown. The prevelent theory 
for the onset of breakdown is principally due to Benjamin* 10 '. In Benjamin’s theory, breakdown is 
explained as a transition between a supercritical upstream flow incapable of supporting upstream 
propagating waves and a subcritical flow which allows upstream and downstream propagating 
waves. The theory is based on an inviscid quasi-cylindrical approximation which neglects all 
axial gradients. In general, no consistent correlation has been found between the occurence of 
breakdown in the numerical solutions of the Navier-Stokes equations and the criticality condition 
of Benjamin. Comprehensive reviews of the progress made in understanding and predicting vortex 
breakdown have been given by Hall* 11 ' and Leibovich* 12 ** 13 '. 

Numerical simulation of the bubble type breakdown of isolated vortex has been simulated by 
many researchers. Grabowski and Berger* 14 ', using Chorin’s artificial compressibility method and 
a primitive variable formulation, calculated the first numerical solution of the steady axisymmetric 
Navier-Stokes equations for this problem. However, the numerical scheme had difficulties in 
realizing fully-converged solutions and obtaining solutions for Reynolds numbers, based on vortex- 
core radius, greater than 200. Hafez et al.* 15 ', using upwind differences and vertical line relaxation, 
solved the Navier-Stokes equations using the streamfunction - vorticity formulation and confirmed 
the earlier results of Grabowski and Berger. Although they were able to obtain converged steady 
state solutions for low Reynolds numbers, they encountered the same difficulty as Grabowski and 
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Berger beyond a Reynolds number of 200. Krause et al. [16J and Menne [17] , used an Alternating 
Direction Implicit (ADI) procedure to solve the same equations as reference 15, but failed to 
obtain a steady state solution. In order to overcome the 200 Reynolds number barrier, Hafez et 
al. tl8J , Salas and Kuruvili [19] and Beren [20] , each working independently, attacked the problem 
using direct matrix inversion technique. In reference 19, the present authors give details of this 
scheme and report solutions obtained for Reynolds numbers as high as 1800. Encouraged by the 
success of this technique, this effort was directed towards finding the effect of a perturbation on the 
spatial stability of these solutions and to eventually capture the spiral-type breakdown mode. The 
full three-dimensional Navier-Stokes equations in a cylindrical coordinate system, decomposed in 
the azimuthal-direction using Fourier transforms, was solved for the first two Fourier modes. In 
this study [21][22] is was found that, at low Reynolds numbers, the effect of perturbation on the first 
non-axisymmetric Fourier component is relatively small. However, the effect was significant at 
higher Reynolds numbers Inclusion of higher modes, to validate these results, using this method, 
leads to a formidably large and expensive problem. Hence the present effort was directed towards 
solving the full three-dimensional problem using an iterative technique. 

In this study the integral form of the unsteady compressible three-dimensional Navier-Stokes 
equation are discretized in space using the finite volume approach and integrated in time using an 
ADI scheme. The in viscid components of the equations are discretized using Roe’s flux-difference 
split upwind scheme and the viscous components are discretized using central difference. Multigrid 
acceleration is used to accelerate the convergence when the problem has a steady state solution. 
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2. Governing Equations 


The fundamental equations of fluid dynamics, the Navier-Stokes equations, that are based on 
the universal laws of mass, energy and momentum, a constitutive law defining the relationship 
between fluid properties and empirical laws stating the dependence of viscosity and thermal 
conductivity with other flow variables, completely describes all flow phenomenon. Assuming that 
there are no external forces and heat sources, the nondimensionalized Navier-Stokes equations in 
a Cartesian coordinate system can be written as 


8Q *(*-*) , 3 ( 6 ~ g ) 

dt dx dy dz 


= 0 


(2-1) 


where Q is the vector of conserved variables, E,F and G are the vectors of Euler (advective) 
fluxes and £, T and Q are the vectors of viscous (diffusive) fluxes. The equations are arranged 
such that the first and second rows of the vector correspond to the continuity and the energy 
equations respectively, while the third, fourth and the fifth rows correspond to the momentum 
equations in the x, y and 2 directions respectively. Thus 
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where p is the density, e is the total energy per unit mass, u, v, w are the velocities in the x , y, z 
directions respectively and p is the pressure . H, the total enthalpy, is given by 

H = e+- (2-6) 

P 

Tij are the components of the shear-stress tensor, q, are the components of the heat-flux vector 
and t is the time. Assuming that the fluid is Newtonian and invoking the Stokes hypothesis, the 
components of the shear- stress tensor are given by 


2 Moo ( 

du 

dv 

dw 

Txx ~ 3 Re^ 

dx 

dy 

dz 

2 Moo ( 

0— 

du 

dw 

Tyy ~ 3 Re H 

~dy 

dx 

dz 

2 Moo ( 

dw_ 

du 

dv 

T " = 3 flcH 

dz 

dx 

dy 


_ M oo f du , jh\ 

Txy-Ty* ~ Re Py dy + dx J 
Moo ( 9 w d u \ 

Txi==Txx = lh fX {d^ + d^) 

Moo (dv dw\ 

Tyz ~ Tzy ~ ~te f \dz + dy ) 

and the components of the heat-flux vector are given by 
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( 2 - 14 ) 
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where p is the coefficient of molecular viscosity, T is the temperature, 7 is the ratio of specific 
heats, Pr is the Prandtl number, Re is the reference Reynolds number and Moo is the reference 
Mach number. 

The quantities used for nondimensionalization of the equations are a reference length L 
and reference values of density poo, speed of sound a 00. and viscosity p^. Using these, the 
nondimensional quantities are 
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where the dimensional quantities are denoted by a hat and (R is the gas constant. The reference 
Reynolds number is given by 


Re = Poo u o o_L 

P OO 


and the reference Mach number is given by 


Moo 


Uoo 

Q>oo 


where Uoo is a reference velocity. The Prandtl number is given by 


( 2 - 17 ) 
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where k is the coefficient of thermal conductivity. The equation set is closed using a constitutive 
relation, the perfect gas equation of state 

p = (7 - \)p e — ^ (u 2 T- v 2 + w 2 ) (2-20) 

The coefficient of molecular viscosity is determined using the power law 


p = T w (2-21) 

where a; is a constant. 

For convenience of discretizing, the governing equations are transformed from the physical 
domain ( x , y, z) to the computational domain (£, 77 , () using the transformation 

t = £(z,y>z) 

r) = T](x,y,z) (2-22) 

C = C(2M/,z) 

Applying the chain rule of partial differentiation, the governing equations in the computational 
domain can be written as 


where 
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(2-23) 
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where 

r y > |Vr| = ( r 2 x + r 2 y + rj) 1/2 (2-25) 

r *. 

where r stands for (£, 77 , (). The metric terms associated with the coordinate transformation are 
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£z — J (yri z £ y<;z v ) 
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where J , the Jacobian of the transformation is given by 

J = i /[ x t(yn z < - y< z v) - x v(yt z c - + x c(yz z v - y^)] ( 2 - 27 ) 

The metrics can be determined either from analytic expressions for the inverse of the transforma- 
tion, 

x = x (t,V>() 
y = y(t,nX) 
z = z (t,*h 0 

or numerically as a direct result of a grid generation scheme, 

x = x(i,j,k) 

y = 

z = z(i,j,k ) 

In this study the metrics are evaluated, numerically, using the coordinates of the grid points. The 
grid generation scheme and the evaluation of the metrics will be discussed later. 


(2-28) 


(2-29) 
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3. Computational Domain and Grid 


In this study, the Navier-Stokes equations are solved, in a space in the shape of a cuboid as 
shown in Fig. 3-1. The boundaries of the computational domain are referred to in the following 
way. In Fig. 3-1, the left hand side face (yz-plane at x = 0) is called the incoming-face and 
the opposite face is called the outgoing-face. The four faces around the x-axis are called the 
side-faces. An isolated vortex is introduced at the incoming face by specifying the profiles of the 
flow variables. The axis of the incoming vortex is aligned along the x-direction and it is centered 
on the incoming-face. The incoming- and outgoing- faces are squares. Note that the boundaries 
around the vortex are permeable. 



Figure 3-1. The Computational Domain 


An algebraic grid that is uniform in the x-direction and stretched in the y- and z-directions is 
used. The grid in the y-direction, for example, is stretched using the equation 


y — ymax@ 


' 1-1 

.P -o 


0<6<l 


(3-1) 


ft is a free parameter that is used to control the amount of stretching. Fig 3-2 shows the grid that 
is used in this study. It has 65x33x33 grid points in a computational space that is 10x8x8. All 
the cells in the flowfield are rectangular boxes. From the numerical grid generation scheme, we 
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define, 


x 
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= n>±i,j 


Z — Z 





(3-2) 


Eq. 3-2 specifies the coordinates of the discrete cells in the flowfield. 



Figure 3-2. The Grid 


For complex topologies a grid generation scheme produces cells that are, at best, hexahedrons. 
Consider a cell (i,j, k ) in the flowfield, shown in Fig. 3-3. Let Q tih k> the conserved variables, be 
located at the cell-center. A semidiscrete finite-volume representation of the governing equations, 
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Eq. 2-23, can be written as 


1 dQ A (E-S) AjF-T) A (G - g) 
J di + At/ AC 


Choosing AC = A^ = AC = 1, the flux balance in a cell can be written as 



+ (E- S) i+ i ] k - (E - 
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(3-3) 


(3-4) 


+(G-Q) itjMr (G-G\ hH = Q 

where the conserved variables, located at the cell-center (z, j, k), are cell averages and the fluxes 
are evaluated at the cell interfaces i ± 5 , j ± 5 and k ± 5 . It can be shown (Appendix- A) that 
Eq. 2-23 is a consistent approximation to the integral form of the governing equations. 



O Location of Coordinates 


• Location of Variables X Cell Interface where fluxes 

are evaluated 


Figure 3-3. A Cell in the Flowfield 
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4. Time Integration 

We are interested in advancing the governing equations in time from some set of initial 
conditions. We are interested in the steady state solution if one exists or in the time evolution 
of the flowfield if the flow is unsteady. We choose an implicit scheme, developed by Beam and 
Warming* 23 ^ to march the governing equations in time. The scheme is implemented as described 
below. The governing equations, Eq. 2-23 at time level n is 


i/ag\ n [ d[E-e) d(F-f) 
J\dt ) \ d£ dr) + 


d(G-g)f 
dC J 


= 0 


(4-1) 


The solution Q n+1 at the next time level, n + 1 can be evaluated using the backward Euler 
formula 


<r +i = <r+Af(|?j 
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(4-2) 


where At is the time increment. Substituting Eq. 4-2 in Eq. 4-1 we get 
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Using Taylor-series expansion we can write 
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Substituting Eq. 4-5 in Eq. 4-3 we get 
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where I is an identity matrix, and 6^ are difference operators in 77 and C directions 

respectively, 


R n 
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(4-7) 


is the steady state residual of the inviscid part of the governing equations and 
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is the steady state residual of the viscous part of the governing equations. If the flow has a steady 
state the total residual R n - 7l n goes to zero as time goes to infinity. The left hand side of Eq. 4-6 
can be approximately factored as follows 
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Eq. 4-9 can be replaced by the following alternating direction sequence. 


1+JAtSr, 

]+JAt6 c 



A Q = - JAt(R n -7l n ) 
A Q = A Q 
A Q - A Q 


(4- 10a) 
(4- 10b) 
(4- 10c) 


Each of these equations is a set of block tri- or penta- diagonal equations depending on the spatial 
accuracy of the left hand side operator. These equations can be solved sequentially to obtain A Q 
from which the solution tit the next time level can be obtained as 


gn+l = gn +A(? (4-11) 

The convergence rate of this scheme is sensitive to the factorization error. At large time 
steps, this error dominates the left hand side and the amplification factor of the error tends to 
unity. However the rate of convergence is relatively insensitive to the accuracy of the spatial 
differencing on either sides of Eq. 4-10 [24) . Hence a first order accurate spatial differencing of 
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the left hand side, that results in a block tri-diagonal system of equations, is a good choice 
since the computational work required is greatly reduced. However, higher order accurate spatial 
differencing is required on the right hand side for the accuracy of the solution. 

Spatial Differencing 

The sequence of equations, Eq. 4-10 are discretized in the following way. Dropping the 
superscript n for convenience and neglecting higher order terms, Eq. 4- 10a at a cell i,j, k can be 
written as 
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(4-13) 


The left hand side of the other two equations, Eqs. 4- 10b and 4- 10c, in the sequence can be 
discretized in a similar way. The right hand side of these equations are obtained from the solution 
of the preceding equations respectively. 

In the absence of external or internal sources, a flow configuration is the result of a balance 
between the inviscid (advective) fluxes, due to the convection of the flow, and the viscous fluxes 
due diffusion of the flow. Viscous fluxes influence the whole space domain and are independent of 
the flow direction while, the inviscid fluxes are dominated by waves, and information is transmitted 
through specific regions of space determined by the wave-propagation direction. Mathematically, 
the unsteady Navier-Stokes equations are parabolic when the flow is viscous dominated and 
hyperbolic if the viscous effects are negligible. This difference in the physical behavior between 
the inviscid and viscous fluxes has to be taken into account by the numerical discretization, in 
order to model the problem properly. In the next two sections we describe, how the flux vectors 
E,F,G,£,F and Q and the flux Jacobian matrices dE/dQ,dF/dQ,dG/dQ,d£/dQ,dF/dQ 
and dQjdQ are evaluated at the cell interfaces. 
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5. Discretization of the Inviscid Fluxes 

In this section we outline the discretization of the inviscid fluxes using the Roe s flux-difference 
split upwind scheme. This scheme can be described using the one-dimensional Euler equations 

WJffl (5-1) 

Jdt dt, 

which is a set of three equations. We seek a solution Q(£,t) that satisfies Eq. 5-1 and the initial 
condition 


Gtf,0) = Qo(0 


(5-2) 


A semidiscrete finite-volume representation of Eq. 5-1 can be written as (Refer Eq. 3-3) 



Once the fluxes E at the cell interfaces i + ^ and i — ^ ‘tre evaluated as a function of Q, Eq. 5-3 
can be integrated in time. 

To determine the fluxes at the cell interface, we use Roe’s flux-difference splitting upwind 
scheme [25] Both analysis and numerical experiments have shown that this scheme provides a 
more accurate representation of shock waves and boundary layers than other conservative upwind 
schemes. However, the computational work required is slightly greater than that of other upwind 
schemes. Roe’s scheme constructs the solution to Eq. 5-3 by solving an approximate Riemann 
problem at each of the cell interfaces. We begin by constructing a higher order interpolation 1261 
of the state variable Q on either sides of a cell interface, namely: 


((?,,), - + i = Qi + ^[(1 + K)8Qf + (1 - K)8Q i ] 

1 = Qi-t-i — 4 [( 1 + K )^»+i + 1 1 ~ K )$Qi+ 1] 

where 

6Qt = Qi+i~Qi 

and 

6Q- = Qi- Qi - 1 


(5-4) 


(5-5) 


(5-6) 
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K 

Scheme 

2nd Order Truncation Error 

-l 

Fully Upwind 


0 

Fromm’s 

d*Q/de) 

1/3 

Third Order 

0 

1 

Centered 

-jAfVWO 


Table 1. Values of k and Associated Truncation Error^ 2 *^ 


The parameter k determines the accuracy of the spatial differencing. The common choices for k 
and their truncation errors are given in Table 1. For all values of k the scheme is second order 
accurate in space except for k = 1/3 when the scheme is third order accurate (strictly valid for a 
one dimensional problem). When higher-order differencing is used, the gradients SQ have to be 
limited in regions of flow with discontinuities, such as shocks, to avoid oscillations of the flow 
variables. The choice of limiter is postponed until later. 

If E is a linear function of Q, i.e., dE/dQ is a constant, then given the initial left and right 
states by Eq. 5-4, the exact solution of the Riemann problem can be written in terms of the flux 
difference as 


Er-E l = K t(Q R - Q l ) 


(5-7) 


where 




dE 

dQ 


(5-8) 


if the flux Jacobian matrix. For the Euler equations it is possible to write the flux Jacobian matrix 
as 1271 


K t = 1 (5-9) 

where columns of are the right eigenvectors of rows of T~ 1 are the left eigenvectors of 
and is a diagonal matrix containing the eigenvalues of K^. The evaluation of T$, T7 1 
and will be discussed later in the section. Substituting Eq. 5-9 in Eq. 5-7 we get s 

E r -E l = T f A (Tj-'Wr - Q t ) (5.10) 
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(5-11) 


Carrying out the matrix multiplication on the right-hand-side of Eq. 5-10, we can write 

3 

Er - = a ^e n 

n — \ 

where a„ is the projection of the difference between the initial left and right states on to the 
left eigenvectors of K^, i.e., o n is the n* element of the vector T ^ 1 {Qr — Ql)> ^ n is the nth 
eigenvalue and e n the n ti right eigenvector of K^. Each term of the summation represents the 
change in flux caused by the corresponding wave. a n represents the strength of the n* wave, and 
A„ is the wave speed. This is shown schematically in Fig. 5-1. From Fig. 5-1, it is evident that 
the flux at the interface can be computed from either 

E x+ i(Q l ,Qr) = E L + ]T> n A„£ n (5-12) 

or 

+ 

E i+ L(Q L ,Q R ) = Er - J> n A n e„ (5 ' 13) 

— + 

where and denote ihe summation of the change in flux due to negative and positive wave 
speeds respectively. Averaging Eq. 5-12 and Eq. 5-13 we get* 25 ! 

if 3- 

i’,- + i(Ql,£?rt) = g e l + E R -^2 a n\K |e« (5-14) 

L n— 1 

Eq. 5-14 can be written in the matrix notation as 

E,^(Ql,Qr) = \[El + E r - |K { |(Q fl - Oi)] (5-15) 

where 

|K C | = T^|A f |T“ 1 (5-16) 

For Euler equations, where E is a nonlinear function of Q, we can define a locally constant 
matrix whose eigenvalues and eigenvectors satisfy Eq. 5-11 and also the condition 

3 

qr- Ql = Y! ( 5-17 ) 

n= 1 
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For this method to return the exact solution whenever the left and right states lie on opposite sides 
of a shock wave or a contact discontinuity, the Rankine-Hugoniot condition, 

E r -E l = S(Qr-Q l ) (5-18) 

where 5 is the speed of the discontinuity must be satisfied. K^, called the Roe-averaged matrix, 
can be written in the form 


K* = T^T' 1 (5-19) 

where is the right eigenvector matrix of K^, T^ 1 is the left eigenvector matrix of and 
n c is a diagonal matrix containing the eigenvalues of K^. Hence the flux difference in Eq. 5-15 
can be written as 

E, H (Qi,Qr) = 1[E L + Ej ) - |K £ |(Qk - Qi)] (5-20) 


where 


K £ | = T { |n £ |T £ -' 


(5-21) 



Figure 5-1. Schematic Representation of Waves at an Cell Interface 
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The matrices T^, T^" 1 and II^ are obtained by replacing the quantities in T^, 1 and by their 

Roe averages. The Roe averaged quantities are evaluated as follows 1281 

P* = \/PLpR (5-22) 


r* - VEMi + 

\fpl + y/pR 


(5-23) 


where 

f = (u,v,w,H ) (5-24) 


Again, the evaluation of T^, 1 and will be presented later in this section. 


Limiter 


There are several choices for the limiter 1291 to ensure that the Total Variation of the initial 
distribution of Q is Diminishing (TVD). Here we choose the minmod limiter where 8Qf and 8Q^ 
are given by the following expressions: 


8Qf = minmod[(Q,- f i - Qi),/3(Qi - Qi-i)] 
8Q~ = minmod[(Q t - Qi-\)-,P{Qi+\ - Qi)} 


(5-25) 


where 


minmod 



0 

sign(x)min( |x|, |y|) 


sign(x ) / sign(y ) 
^iyn(:r) = sign(y ) 


(5-26) 


and 


0 = 


3 — K 

1 — K 


(5-27) 


is a compression parameter. In smooth regions of the flow, the effect of the limiter is negligible. 
However, in the regions of the flow with high gradients, the accuracy reduces to first order. 
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Three-dimensional Problem 

The one-dimensional flux differencing technique can be extended to multidimensions by 
treating each direction in a locally one-dimensional manner. Although this approach allow 
interactions, of the flowfield, to occur only normal to the cell faces, it has been very successful in 
solving a variety of problems. For the three-dimensional problem the interface fluxes in each of 
the three directions are written as 

E i+li,k = 2 [^<0 + E R(i) - I K d (<?*(•■) “ #£(•■))] 

= 2 ^ 0 ) + Fr u) ~ ~ Qm)) 1 < 5 ' 28 ) 

G i,j,k+l = 2 [^(*) + G R(k) - |K c |(Q7t ( *) - QL(k))] 

Here the flux-vectors have 5 elements and the flux Jacobian matrices are 5x5. 

The contribution of the inviscid terms to the left hand side of Eq. 4-13 are evaluated as follows. 
The flux at the interface i + k is 

E i+$,j,k = \[ e L + Er~ - Q l )\ (5-29) 


Linearizing the flux at this cell-interface, we get 


d Qi+kJ,k 

Carrying out the differentiation we get 


dE i+ l - i k dE i+'- i k d£,+i ? k 

v + ^^-AQ R 


dQi 


OQr 


(5-30) 


(™aq) = I 


9 El A Q l + ^AQr-Ik^AQr-AQl) - ^A Q(Qr-Ql) 


OQl 


dQ R 


dQ 


(5-31) 


Neglecting the higher-order term (d|Kf|/d<2) AQ(Qr - Q L ) we can write 

(f A<? ), +u * k H(f a< 3 X + (l A0 )* - W(A<3 * - A<w 

Making the approximation ( dE/dQ) L = ( dE/dQ) R = we can write 

'3£) \ K, + |K,I._ K,-|K 


( 


K*+ K f Kc 

2 AOt + — T 


AQr 




(5-32) 


(5-33) 
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where IT, and 11^ are diagonal matrices containing the positive and negative eigenvalues respec- 
tively. Note that the matrices are evaluated at the cell interface using the Roe-averaged quantities. 
In order to obtain a tridiagonal system of equations, we choose Ql and Qr to be equal to Qi,j,k 
and Qi+i t j t k respectively, which is a first order approximation. Thus Eq. 5-33 becomes 




(5-34) 


Similarly we can write 




(5-35) 


Eigenvalues and Eigenvectors 

Here we describe the evaluation of the eigenvalues and eigenvectors of the inviscid flux 
Jacobian matrices dE/dQ, dF/dQ and dG/dQ. The inviscid fluxes E, F and G can be written 
in a general form as 


where 


E,FovG = 


| Vr 

J 


pu 

pHu 

puu + pr x 
pvu + pf y 
pwu + pf z 


for r = £,rj or £ 


(5-36) 


u = ur x + vr y -f wr z (5-37) 

is the contravariant veloci ty component in the r-direction where r stands for £, r/ or (. Let 

K( = OE/dQ 

K, = dF/dQ (5-38) 

K< = dG/dQ 

The inviscid flux Jacobians K^, K,, and can be written in a general form as 


0 

0 

r x 


r z 


(aq^ — H)u 

7 XL 

Hr x ~cruu 

Hr y — avu 

Hr z —crwu 


aq 2 r x —ui 

of x 

uf x —&ur x + u 

uf y — avr x 

uf z — awr x 

(5-39) 

aq^fy — vh 

Ofy 

vr x — our y 

vr y —avf y -\-u 

vr z —uwfy 


_ a q i r z —wi 

of z 

wf x — auf z 

wf y — avr z 

wfz—owfz+u _ 
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where 


9 2 = 


2 i 2 I 2 

U + V + W 


( 5 - 40 ) 


a = 7 — 1 


( 5 - 41 ) 


and r stands for rj or (. It is difficult to determine the eigenvalues of K r . However by similarity 
transformation we can obtain a matrix P r whose eigenvalues are easy to evaluate. We can write 


P r = M _1 K r M 


( 5 - 42 ) 


where 


M = 


dQ 

dU 


Q 

u 

V 

w 


0 0 0 0 
\j a pu pv pw 
0/9 0 


0 

0 


0 

0 


p 

0 


0 

0 

p 


( 5 - 43 ) 


and the inverse 



■ 1 

0 

0 

0 

0 ‘ 

Vk 

Jl 

7 

aq 2 

a 

-au 

— (TV 

—aw 

-u/p 

0 

l/P 

0 

0 

-v/p 

0 

0 

l /p 

0 


.-V>/P 

0 

0 

0 

!//>. 


( 5 - 44 ) 


and 


U = 


P 

P 

u 

V 

w 


( 5 - 45 ) 
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Using Eqs. 5-39, 5-43 and 5-44 we get 



0 pr x 

u pa?r x 
f x lp u 

r y /p 0 

r z /p 0 


pr y pr z 

2 - 2 - 

pa L r y pa*r z 

0 0 

u 0 

0 u 


(5-46) 


Since matrices K r and I* r are similar, their eigenvalues are the same^. Hence the eigenvalues 
of K r can be found by determining the eigenvalues of P r . The eigenvalues of P r are 


1 

-ft 
1 


m u + a " 

a; 


u — a 

A°r 

_ |Vr| 

u 

J 

A° r 


u 

1 


u 


A+ and A r are referred to as the acoustic eigenvalues and is referred to as the convective 
eigenvalue. 

Now we seek a matrix S r that will diagonalize P r such that 

S r _1 P r S r = A r (5-48) 


where 


A r 


'A+ 0 0 0 O' 

o A; o o o 

0 0 A° r 0 0 

0 0 0 A° r 0 

. 0 0 0 0 A° r _ 


(5-49) 


Columns of S r are the right eigenvectors of P r and the rows of S7 1 are the left eigenvectors of P r . 
Skipping the algebra, an eigenvector matrix and its inverse that will diagonalize P r as prescribed 
by Eq. 5-48 are 
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(5-50) 


and 



' 1 

1 

1 

0 

0 ‘ 


a 2 

a 2 

0 

0 

0 

s r = 

ar x /p 

~ar x /p 

0 

h/p 

fh x f p 


ar y /p 

-af y jp 

0 

[y/p 

ffly/p 


. af z/p 

~ar z /p 

0 

hip 

fh z /p 


'0 

l/2a 2 

pr x / 2a 

pr y j 2a 

pr z /2a ' 

0 

l/2a 2 

-pr x /2a 

-pr y /2a 

—pf z /2a 

1 

— 1/a 2 

0 

0 

0 

0 

0 

ph 

Ph 

Ph 

0 

0 

pfh x 

pfhy 

pfh z 


(5-51) 


where l x , l y , l z and fh x ,m y ,m z are the components of two, arbitrary, mutually perpendicular unit 
vectors I and m respectively that are orthogonal to r whose components are f x ,f y ,f z . 

Using Eqs. 5-42 and 5-48 we can write 


K r = T r A r T r - 


-l 


(5-52) 


where 


and 


T r = MS r = 


1 


1 

1 

0 

0 “ 

H + au 

H 

— au 

q 2 

V 

w 

u + ar x 

u - 

- af x 

u 

/« 

fh x 

v + af y 

v - 

-afy 

V 

h 

Thy 

w + af z 

w - 

- af z 

w 

iz 

fh z _ 


(5-53) 




where 


1 

2a* 


and 


aq* — au 

a 

ar x — au 

af y — av 

af z — aw 


0 _ 
aq +au 

a 

—ar x — au 

—i afy — av 

—af z — aw 


2a 2 — 2aq* 

—2a 

2au 

2av 

2 aw 

(5-54) 

— 2a 2 i> 

0 

2a* l x 

2(1* ly 

2a* l z 


—2a*w 

0 

2a*fh x 

2a*fh y 

2a*fh z 



V 

= ul x + vl y 

+ wl z 


(5-55) 


w = 

uih x + vfhy 

+ wfh z 


(5-56) 
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are two mutually perpendicular velocity components that are orthogonal to u. Substituting Eqs. 5- 
49, 5-53 and 5-54 in Eq. 5-52, we can obtain the matrix K r . Matrix |K r | can be obtained from K r 
by replacing all the flow quantities by their Roe averages given by Eqs. 5-22, 5-23 and 5-24, and 
replacing the eigenvalues by their magnitude. The flux differences |K r \{Q r — Ql) in Eq. 5-28 
can be written in a computationally efficient form as 


|K r |(Qfi-Oi)=^|u*+a , | 

J 




H* + a*u* 
u* + r x a* 
v* + r y a * 
w* + f z a* 

1 

H* - a*u* 
u* — f x a* 
v * — f y a* 
w* — f z a* 


I 


* L >' 


1 

m 2 


(«*) 


u 


w 


PR ~ PL 


PR ~ PL 

P*(UR-U L ) 

i(a-f 

+ 2a* 

PR ~ PL 

p*(u R - u L ) 

2(a*) 2 

2a* 

PR- PL 
(a*) 2 

} 


|Vr| 


«i*l P' 


U*(UR -u l) + V*(vr - Vl) -I- W*(WR - Wl) - U*(UR ~ Ur) 
UR - U L - r x (y,R - Ui) 

vr-vl- ry(-UR - ul) 

wr - wl - r z (uR ~ ui) 


(5-57) 


where the quantities with a star are Roe-averaged quantities. Using Eq. 5-57 in Eq. 5-28 the 
contribution to the right hand side of Eq. 4-13, due to the inviscid terms can be evaluated. In 
order to determine the left hand side, we have to evaluate the matrices T r n^ T r and T r II r T r 
The evaluation the matrices Trll+T^ 1 and TrlljrT;: 1 is given in Appendix B. 
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6. Discretization of the Viscous Fluxes 

In this chapter we discuss the treatment of the viscous terms. The viscous fluxes £, T and Q 
can be written in a general form as 


£, T or Q 


jVr[ 

J 


0 

UT r { + VT r1 , + WT r ( - q T 

T r( 

Tr V 

TrC 


for r— £,t] or £ 


( 6 - 1 ) 


where 


T t £ = Txx'f'x "I- Tyx^y H“ T zx f z 

T ri) — T X yr x + Tyyfy + T zy^z ( 6 - 2 ) 

T r { = T xt r x -f- Tyzfy + Tzz’l'z 


and 

9r = q x r x + q v f y + qzf z (6-3) 

where r stands for £, r) or (. Note that the shear-stress and heat-transfer terms contain spatial 
derivatives with respect to x, y and z. They can be evaluated in the transformed space, (£,/?, C)> 
using the chainrule and hence we can write 


u x 

V x 

W x 

Tx 


£x 

Vx 

Cx 



v t 

W£ 

T 0 

Uy 

Vy 

Wy 

Ty 

= 


Vy 

Cv 



V„ 

W v 

Ty 

. u * 

V z 

Wz 

Tz\ 


u* 

Vz 

C*J 


.«C 

V C 

Wq 

*cJ 


Substituting for the shear-stress and heat-transfer terms in Eq. 6-1 from Eqs. 2-7 through 2-15 and 
using Eq. 6-4 we can write the viscous fluxes as 


where 


£, T or Q = 


r o i 


|Vr| Mqq 
J Re ^ 


$ 2 r 

* 4 r 


* 5 r 


for r— (,ij or £ 


(6-5) 


$3r = <t>l( + + 4> + <f> 2yrVy + fayrMy + ^lCr^C + <f>2(rV<; + $3Cr w < 

$ 4r = <f>4{ r U{ + <j>ft£ r V£ + <j)^ T W£ + (f>4 vr U, + <f> ^ r V v + W n + + <j> 5( r V( + 

$f)r = + 4>&(rV( + + <f>7yr^y + ^ T V n + fa^Wy + <j> 7 ( r ll( + <f>^ T V( + <t>c% T W(; ^-6) 

^ 2 r = + $ 4 *-^ + jrTf) + 
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where 


^ 0 £r = £x^x + £y^y + z 
^Oyr = "Hx^x ~^~ "h Vz^z 
^0 (r ~ (xf x “f“ Cy^y “t“ CzG 


^lfr — * ^0£r + g£x r x 

^ 2 £r ~ £x?y — -jiy^x 
2 

^ 3 £r ~ £x r * — gG r x 

4^\t}t — ^0 rjr + gf/x^x 

, _ - 2 _ 

<P2i?r — ^x^y — ^ r !y r x 

$ 3 fjr — 'Hx'f'z ~ <^Hz r x 

~ + gCx^x 

2 

^ 2 Cr = Cx^y ~ gCy^x 
2 

^ 3 Cr = Cx r * “ gG r x 


^4(r — £y r * ^ xV y 
fatr = <f>i) {r + g^y 

^6f r ~ Cy 7 ^ — gG^y 

, _ - 2 - 
<PAr}T — r 1y r x — g^x^y 

^Stjr = ^Or/r + g* 7 y^y 

<?W = * 7 yG - ~rj z fy 

^ 4 Cr = Cy^x — g Cx^y 
^ 5 (r — ^OCr + gCy^y 
^f>C r = Cy^x “ g G ^y 


(6-7) 


( 6 - 8 ) 


(6-9) 
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^7£r — £z r z z 

= tzfy ~ ^yfz 

<t>^r = 0O£r + 3^^ 

, _ - 2 _ 
<P7r)r — Vz r x <^)xT z 

hvr = Wy - | r?j,r z 

*^9t jr = ^Ojjr "H qVz^z 

<^7£r — Cz^x — gCx^z 

^8Cr = Cz^» — g Cyi'z 

4>^r = </>0Cr + gC*r 2 


( 6 - 10 ) 


where r stands for £, r) or (. 

Once the metrics Eq. 2-26 are determined at each cell-interface (Evaluation of metrics will 
be discussed in the next section.), the quantities given by Eqs. 6-7, 6-8, 6-9 and 6-10 can be 
determined and hence the viscous fluxes through the cell-interfaces can be computed from Eqs. 6- 
5 and 6-6. For example, 


(6-11) 

where 

+ HtfWt + <t>h)tUri + + 4>^W V + (f>l^U^ + <f>2(£V( + <f) 

$4£ = 4> 4(t u t + <t> 5(( v £ + + <t>Ai)i u v + <t>hriiV v + + <^4« U C + ^f>Cf 

= ^ 7 ««£ + <f> 8((V ( + fa ((W( + fa rii^n + <f> * v t v v + <l%t w n + fact u C + fa(t v C + fau w <i ^ 1 ^ 

+ $4£^ + + ~^p + fa v {T v + fa^U 



where 
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U t 

V ( 

T* 


~ t - fi+l,j,k fi,j,k 


(6-13) 




u r 


tU* 


— (fi])i + L t j t k ~ Afi+l,j+\,k + fi,j+l,k fi+l,j-l,k fi,j-l,k) 


(6-14) 


L ^*7 




U C 

V C 

T C . 


” 4 (/•+!, i,*+l + /tj.fc+l /«+!,>,*- 1 




and 


(6-15) 


u 

V 

U 7 





i+$,j>k 


2 (f*+hj,k + /«,>,*) 


(6-16) 


Notice that the cross-derivative terms are evaluated using a symmetric difference molecule and 
hence there is no contribution from these terms to the diagonal element of the coefficient matrix. 
Similarly the viscous fluxes at the other cell-interfaces can be determined. From the viscous fluxes 
at the cell-interface, the contribution of the viscous terms to the right hand side of Eq. 4-13 can 
be determined. 

On the left hand side we have to evaluate the flux Jacobian matrices d£fdQ,dTjdQ and 
dQ jdQ of the viscous fluxes S, T and Q. In order to obtain a tridiagonal system of equations, 
we neglect the cross-derivatives before the fluxes are linearized. For example, consider the flux- 
vector £. First we neglect till the derivatives with respect to rj and (. Then we discretize the terms 
according to Eqs. 6-13 and 6-16 and linearize them in the following manner. We can write, 


d£ a ~ 


dS i+ i jk d£ i+ i k 

= J,k + 

dQi+l,j,k d Qi,j,k 


(6-17) 


Assuming that the molecular viscosity n is locally constant and carrying out the differentiation 
and simplifying we can write 


d£ . \ 

a Q A + j J,„ 


( 6 - 18 ) 
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where 


and 


frWffiL, 

(v-) =-(^] 

' *'*+| j .* \ d Q/i,j,k 


(6-19) 


( 6 - 20 ) 


The viscous flux Jacobian matrices ( Vt ) and f V, ) are functionally the same. The 

flow variables in (V^J ] except n are evaluated at the (i + l,j, k) ih cell-center and the flow 
variables in ( V7 ) except /i are evaluated at the (i,j, k) th cell-center. However the metrics 


are evaluated at the (i + ^,j, k) ih cell-interface. \i is evaluated at the cell-interface as an average 
given by Eq. 6-16 Similarly, 





( 6 - 21 ) 


where 


and 



( 6 - 22 ) 



(6-23) 


The other viscous flux Jacobian matrices dT /dQ and dQ / dQ can be evaluated in a similar manner. 
Let V r denote the general viscous flux Jacobian matrix where, r stands for £,77 and C such that 


V* = dt/dQ 

V, = dT/dQ (6-24) 

V c = dQ/dQ 


V r can be written in a general form as 
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X 


V r = J 



Moo A 

Re (.• 


0 

-2g 2 — 3u 2 — 7 (e — 2g 2 ) jPr 

—u — f x u 

— V—fyU 

—w — f z u 


0 0 

7 jPr u + f x u— , yu/Pr 
0 1 + f x r x /3 

0 rj,f x /3 

0 r z f x / 3 


0 

v + fyU — 'yv/Pr 

Tx 1'y /3 
1 + r y r w /3 
r z f y /3 


0 

w + r z u — ~fw/Pr 
r x fz/ 3 

f y f z /Z 

1 + ^r*/3 


(6-25) 


where 


u = -^(ur x + vr y + wr z ) 


(6-26) 


Again all the flow variables except /i are evaluated at the cell-center, fi is evaluated as the average 
of the adjacent cell values (Eq. 6-16). All the metrics, area and volume, 1/J, are evaluated at the 
cell-interface. The volume at the cell-interface is evaluated as the average of the volume of the 
adjacent cells. Some of the details of the derivation of V r arc given in Appendix C. 
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7. Metrics, Areas and Volume 


From the numerical grid generation scheme, we get (Eq. 3-1), 

* = *(i±i,;±i *±!) 

= (7-1) 

2 = z ( i ± — , ? ± — , A: ± — 

V 2 2’ 2 

Eq. 7-1 describes a set a of hexahedrons in the flowfield. Note that in general, a grid in an arbitrary 
space may not yield hexahedrons, i.e., four adjacent points may not be coplanar. 

Let ABCDEFGH shown in Fig. 7-1 represent a cell (i,j, k) in the flowfield. The cell average of 
the flow variables, Q are located at the cell center (i,j, k). The vertices of the cell are 



Figure 7-1. A Cell 
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A : 
B : 
C: 
D : 
E: 
F : 
G : 
H : 


(•' + 5 >> + 5 ,* + i) 

(i + j.i-j.i + i) 

(> - |.J + 1*- 1) 

(•-j.i + j.t + i) 

0 _ l’ J - 5’* + 1) 


(7-2) 


We need to evaluate the metrics, cell-face areas, cell-face normals and volume, of each of these 
cells, in order to solve the governing equations. Consider the constant £ face ABCD. Let a and b 
be the two diagonal vectors on the face. Then it can be shown that 


kr + k I+ k^i Sxi 


(7-3) 


is the directed face area and the resulting vector is perpendicular to the face. Therefore the unit- 
normal to the face is 


n 


= & + £yj +£zk = 


axb 

axil 


The magnitude of the face area is 


(7-4) 


J 



(7-5) 


Similar evaluations can be done for the other faces. To evaluate the inviscid fluxes on a constant 
£ face, the metrics associated with the other directions namely,^, i/y, i \ z , £ x , £ y and £ z are not 
need at the face. Similarly for the other directions. However all the metric quantities are needed 
at any face in order to determine the viscous fluxes. Evaluations of those metrics that cannot be 
computed using the above, formulas will be discussed later. 
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Figure 7-2. Cell Volume Decomposition 

The volume of the hexahedron is evaluated as follows 1311 . The hexahedron can be decomposed 


into three pyramids as shown in Fig. 7-2. The bases of the three pyramids are ABCD, ABFE and 
ADHE and they have a common apex at G. The volume of each pyramid is one-third the height 
times the area of the base. Therefore the volume of the hexahedron can be written as 


1 

J 




(7-6) 


The notation AG, for example, refers to the vector from A to G pointing towards G. The vectors 
on the pyramid bases are defined such that the directed face areas arc pointing into the cell. This 
ensures that the resulting volume is a positive number. Since a hexahedron has four diagonals, 
its volume can be evaluated four different ways depending on the choice of the diagonal. In this 
study the average of these four volumes is used as the volume of the hexahedron. 

Now we present the evaluation of metrics that are needed to compute the viscous fluxes but, 
cannot be computed using the formulas given above. Let us consider the constant £ face i + ^ , i, k . 
The metric quantities r\ x j J, r] y /J, r\ z /J , C*/ 7, ( y /J and ( 2 / J at this face are evaluated as an 
average of the neighboring face values. Thus 
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where 8 stand for x,y or z. It is also necessary to compute a volume associated with the cell 
face in order to evaluate the viscous terms on the left hand side of the governing equations. This 
volume is evaluated as an average of the neighboring volumes. Thus 


t Am 

'•+5 fi >* — 2 


where J is the reciprocal of the volume. Similarly the metrics on the other cell-interfaces are 
evaluated. 
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8. Boundary Conditions 

A major hurdle in simulating vortex breakdown, especially the spiral-type breakdown, is 
the numerical modelling of boundary conditions consistent with experiments. In the past we 
were able to simulate bubble-type breakdown by specifying that the flow is essentially inviscid 
near the outer edge of the computational domain 1221 . Other investigators have simulated similar 
bubble-type breakdown by specifying pressure distribution along the outer edge 1321 . In the present 
study, we assume that the outer edges of the computational domain are reasonably far away 
from the breakdown region. Hence the flow in this region is considered to be inviscid, steady 
and reasonably benign. The boundary conditions that are used are based on the characteristic 
variables 1331 . These boundary conditions are developed by assuming that the flow is locally 
one-dimensional and perpendicular to the boundary, i.e., the derivatives along the surface of the 
boundary are neglected. With the above assumptions, the following characteristic form of the 
Euler equations can be derived at a constant r boundary, where r stands for f , r/ or (. 


ds 

dt =0 

along 

dr 

— — u 
dt 

(8-la) 

^=0 

dt 

along 

dr 

di = U 

(8- lb) 

— ± - o 

dt pa dt 

along 

(dr Y* 

( — = u ± a 

\dt ) 

(8-lc) 



Figure 8-1. Characteristic Waves in a Flow 
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where s is the entropy, l and V are the component of velocity normal and tangential to the 
boundary respectively. The normal velocity is (Eq. 5-36) 

u = ur x + vv y -f wf z (8-2) 

A sketch of the characteristic paths is given Fig. 8-1. According to Eq. 8- la and Eq. 8-lb , the 
entropy and tangential velocity, respectively, are convected along the direction of u. Assuming 
that the flow is locally homentropic, Eq. 8-lc. can be written as 

dR ± n , /drf 

7t = 0 along (*)=«*“ < 8 - 1 2 3 4 > 

where 

mi - I 2a 

R* = u ± — (8.4) 

a 7 

are the Riemann variables. The boundary conditions are applied by specifying appropriate 
conditions at a set of phantom or ghost cell around the boundary. The characteristic equations 

are used to update the variables at these ghost cells. The ghost cells are updated explicitly at the 

end of each time step. 

For this specific problem the incoming-face is either subsonic or supersonic inflow and the 
outgoing-face is either subsonic or supersonic outflow. Conditions at a side-faces could be any 
one or a combination of the above. Therefore the boundary conditions described are specialized 
for the following four cases. 

1. Supersonic inflow 

2. Supersonic outflow 

3. Subsonic inflow 

4. Subsonic outflow 


Supersonic Inflow 

For supersonic flows all the eigenvalues have the same sign. Hence, since the flow is coming 
into the computational domain all the flow variables in the ghost cells are specified. 


Supersonic Outflow 

Again, since the flow is supersonic, all the eigenvalues have the same sign. Since the flow 
is going out of the computational domain, all the flow variables are extrapolated from inside the 
computational domain to the ghost cells. 
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Subsonic Inflow 


In subsonic flows four of the eigenvalues have the same sign and one has the opposite sign. 
For the subsonic inflow case four of the eigenvalues, u + a, u, u and u are positive and one, 
u - a is negative (Eq. 5-46). Fig. 8-2 shows a schematic representation of the characteristic waves 
at a constant r boundary at time level n + 1. At this time level the R + characteristic can be 
determined from the conditions at the ghost cell at time level n and the R characteristic can be 
determined from the conditions at the first interior cell at the same time level. Thus 


R+ = + 

R" = uf - 


2o£ 

a 

2a” 

a 


(8-5) 


where G denotes the ghost cell and I the first interior cell. The two Riemann variables can be 
added and subtracted to determine a local normal velocity and speed of sound at the time level 
n -j_ This new normal velocity and the speed of sound is assigned as the ghost cell value at time 

level n + 1 . Therefore we can write 

-n+l _ R+ + R ~ _ gfiiJg + “o - (8-6) 

“o - 2 ~ 2 a 


t 



Figure 8-2. Characteristic Waves at an Inflow Boundary 
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(8-7) 


!»-+! - 


<r(R+ - R") _ a(ufi - Sf) a G + a? 


+ 


2 


Since the tangential velocity and entropy are convected in the direction of u, they are held fixed. 
Therefore 


V G +1 = Vfi 


( 8 - 8 ) 


and 


„n-|-l _ n 
S G - Sq 


(8-9) 


As shown in Appendix-B 


vl x + wfh x — u — uf x 

vly + Wfll y = V — ufy 

vl z + wfh z = w — ur z 


( 8 - 10 ) 


are the components of the tangential velocity V in the x,y and 2 directions. Therefore from 
Eq. 8-8 we can write, 


Eq. 8-11 can be rearranged as 


(u - ur *)£ +1 = (u - ur x )n 

(V - Ufy ) G +1 = (V - Ufy) G 

(w - ur*) G +1 = {w - uf z ) r Q 


( 8 - 11 ) 


«a +1 * + («G +1 - 

+ («S +1 - u n a)f y 
«»G +1 = < + (^G +1 - So)r, 


( 8 - 12 ) 


From Eq. 8-9 we can write 


P_ 

pVg 


n+1 / \n 

1 ■ (A 


(8-13) 


From Eqs. 8-7 and 8-13 the total energy per unit mass e and the density p can be easily determined. 
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Subsonic Outflow 

Subsonic outflow conditions can be derived in a similar manner as subsonic inflow conditions. 
Fig. 8-3 shows the characteristic waves at a subsonic outflow boundary. At time level n + 1 the 
R + characteristic is determined from the conditions at the interior cell adjacent to the boundary 
cell at time level n and the R“ characteristic is determined from the conditions at the ghost cell 
at the same time level. Thus we can write 


and therefore 


R + = u? + 
R" = u n G - 


2 a 


» 

I 


a 


2a 


n 

G 


a 


U G ~ 


R- + R+ ufc + - af 


„n+l _ 
°G “ 


Or(R~ — R + ) 


g ( fi G~ fi f) , «G+ q ? 

A ' 


(8-14) 


(8-15) 

(8-16) 


The tangential velocity and entropy are extrapolated to the ghost cell, using zeroth order extrapo- 
lation, from the interior cell. Therefore 


t 



Figure 8-3. Characteristic Waves at an Outflow Boundary 
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( 8 - 17 ) 


( 8 - 18 ) 


Similar to the subsonic inflow case it can be shown that 


■ V 1 = “” + (“« +1 - B ") 


— v \ + ( U G + ~ u \) r y 

= u;" + (uq +1 — «f)r. 


P'M 


Now we have all the conciitions needed to write the complete set of algebraic equations. 
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9. Solution Procedure 

In this section we describe the procedure for solving the algebraic equations obtained by 
discretizing the governing equations. By substituting for the flux Jacobians from Eqs. 5-33, 5-34, 
6-18 and 6-21 in Eq. 4-13 we get, for a cell ( i,j,k ), 

+ JA({ (T<n+ V) |+| . AQ iJtt + (T«n { -T { ->) . +JJ A<?, +ul 

-( T < n r - ( T < n ?V) H , 


JAt{E i+ i j k F i j_i ik + G i j k+ i G i - k _ i 

~ £ i+i,J,k + S i-lj,k - fx,]+lk + Fij-l,k - GijM'j + g i,j,k-0 

Eq. 9-1 can be written as 

A.{AQi-ij k + B,ZiQ, D, 


(9-1) 


(9-2) 


where 


B i ,i + ;A,((T ! n*T ( -‘) i+{ji -(T ! n ( -T-‘) Hj 

+ ( v t) h ^ + ( v ?)h** ) 

c.= ^{(T^T,-*)^ -(vf) H } 


(9-3) 


are 5x5 matrices and 


D« - - JA H E i +i,j,k - E x-\, h k + F i,m,k - F i,j-i,k + - G i jk _i 

- £ '+\,i>k + £ i -r 2 ,],k — F i,j+k,k + F i,j-k,k — Gi,j,k+\ + &i,j,k-k) 


(9-4) 
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is a vector with 5 elements. By ordering these equations in the increasing 
i direction, we get a block tri-diagonal system of equations that can be written as 

B 1 AQ 1 ,j,k + Ci A Q‘2,j,k = Di for i ~ 1 

Ai AQ t _i + Cj AQi+iyjt = D, for i = 2,7—1 (9-5) 

A/A<2/_i Jik +BiAQ lijik =Dj for i = I 

where i = 1 is the first interior cell and i = I is the last interior cell in the z-direction at any 
j, k location. Similar system of equations can be constructed for the other two directions, j and 
k Eqs. 3- 10b and 3- 10c, of the alternating direction sequence where, the resulting equations are 
arranged in their respective direction. 


Thomas Algorithm 


Each of the block tri-diagonal system is solved using the Thomas algorithm where the solution 
is obtained in two stages In the first stage a set of auxiliary quantities E and F are computed 
using a set of forward recursive formula starting at i = 1. In the second stage the required solution 
is obtained from a backward recursive relation starting at i = I. Skipping the derivation, these 
formula are, for the forward sweep, 


Ei ==-B]‘ 1 C 1 

Fi= : B^Di 


E l :=-(A t E,_ 1 +B<r 1 C, 

F, - (A,E,-_i+B,-) -1 (D;— AiFi_i) 


} for i= 2 , 1 


(9-6) 


and for the backward sweep 

&QL), k - F; 

&Qi,j,k = EiAQj +ljijt H- F t for z= 7-1,1 


(9-7) 


The inversion of the 5x5 blocks in Eq. 9-6 is done by Gaussian elimination. This procedure is 
applied at all the j, k locations to obtain the AQij,k f° r the whole flowfield. 

The solution vector A Qij t k arranged in the j direction at each k,i becomes the right hand 
side of the second sequence and the resulting equations are solved in a similar manner as described 
above. This procedure is repeated for the third sequence. From the final A Q the solution at the 
next step is obtained form Eq. 3-11 which is 


Qij,l ~ Qh,k + A Qi,i,k 


(9-8) 
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For steady state problems, this process is repeated until the desired convergence criteria is satisfied. 
For unsteady problems this process is stopped at the desired time level. 


Time Step 

The time step used to advance the solution is determined as follows. The time step limit at a 
cell is given by the following equation 


where 


i _ , i 


At{ 

1 


A!, “ Jy '*2 

1 _ , 1 
— J hi,k o 


At 


1 1 / 1 1 1 \ 

At l>J>k ~ CFL\Att + At v + At J 

. (|e,+ +(|,|+ 

(W + a) ‘.H.‘( J T ! } . +i4 + (I s ' + 

m + a) -H (^L, + asi + a) -^ (TL, 


and the CFL used in this study is 


(9-9) 


(9-10) 


1 < CFL < 10 


(9-11) 


The flow variables at cell-interfaces are determined by averaging the cell-average quantities from 
the adjacent cells. For steady state problems the solution at each cell is advanced using the local 
time step. For unsteady problems a global minimum time step is used. 

Multigrid Acceleration 

For steady state problems, the rate of convergence of the steady-state residual R - 71 can 
be improved using a full approximation multigrid acceleration technique 134 ^ 35 '. Although this 
acceleration technique was originally developed for elliptic problems, it has been successively 
used for a variety of problems including transonic flow problems that contains shocks and other 
discontinuities. In this technique, starting with the finest grid, the problem is solved on a series 
of successively coarser grids and the correction obtained on each of the coarser grids is passed 
back successively and eventually to the finest grid. The coarse grids reduce the low frequencies 
(with respect to the fine grid) errors quicker than the fine grid. Although, for multigrid method, 
computational cost per iteration is higher than that of single grid methods, the gain in the rate 
of convergence far offsets the additional computational overhead. The implementation of the 
muldgrid acceleration technique is presented below. 
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Let 


£h{Qh) = fh 


(9-12) 


represent the problem that we intend to solve where £/, is the nonlinear discrete operator that we 
obtain by discretizing the Navier-Stokes equations, as described in the previous sections and fh 
is the forcing function associated with the equations. For the Navier-Stokes equations, considered 
here, fh = 0. However we will retain it for generality. The subscript h denotes the finest grid and 
2 h, 4 h etc. denotes successively coarser grids where each coarse grid is obtained by removing 
every other grid point from the preceding fine grid. We seek a solution Q h such that the residual 
fh~ £h{Qh) goes to zero after a number of iterations of the solution scheme (ADI scheme in this 
particular work). After jV A iterations of the solution scheme, on grid h, let Q' h be the solution. 
The error in the solution can be written as 


n = Qh- Q' h (9-13) 

or we can write 

Qh=Q' h +i'h (9-14) 


where uh is the correction that has to be added to the fine grid solution to obtain the exact solution. 
In order to obtain this correction, we solve an equivalent problem on the next coarse grid 2 h. The 
equivalent problem on grid 2 h can be written as 

C 2 h{Q2h) =hh (9-15) 


where, 

hh = l[h - C h (Q' h )] + C 2h ( IQ'h) (9-16) 

2 h \2h J 

h 

| represents an operator hat restricts quantities from grid h on to the next courser grid 2 h. The 
2 h 

exact form of the restriction operator will be presented later. The forcing function f 2h ensures 
that the coarse grid problem, Eq. 9-15, is driven by the residual of the finest grid. Starting with 

1 Q' h as an initial guess, the solution scheme is applied to this new problem and let Q'L be the 

2 h 

solution after N 2 h iterations. Again we can write 

Q'lh — Q'ih + u ih (9-17) 
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and we seek to obtain the correction u^h by solving an equivalent problem on the next coarse grid 
4 h. We construct the equivalent problem 

£«(<&) = /« ( 9 - 18 ) 


on the next coarse grid 4/i where, 

2 h /2 h \ 

A»-£j*(Oa) + £» IQ 2 J (9-19) 

4ft \Ah / 

Several iterations of the solution can be performed to Eq. 9-18. Let Q^' k be the solution after 
Nm iterations. This process can be continued on as many successively coarser grids as desired or 
possible. At the end of this process, we proceed to pass (prolongate) the corrections successively 
from the coarse grid to the next finer grid, using a suitable interpolation function, all the way to 
the finest grid. This is done in the following manner. Let us assume that the coarsest level that 

we down to is 4 h and that we were able to solve the problem exactly, i.e., we found Q'[ h . Now 

2ft 

we prolongate the correction Q'[ h - [Q'^h fr° m this grid to the next finer grid 2 h and the corrected 

4 h 

solution Q'^ on grid 2 h can be written as 

Qlh = Q'ih + *2ft (9-20) 


where 

2h 

V2h = T 

Ah 


where f is a prolongation operator. The functional form of this operator will be presented later. 

Ah 

Starting with the corrected solution Eq. 9-15 can be iterated several times to dampen the errors 
introduced by the prolongation operator and obtain a new solution on grid 2 h. Let Q', lh be the 

h 

solution after A/ 2 h iterations of the solution scheme. The correction Q‘ 2h - [Q' h is prolongated to 

2 ft 

the finest grid and the corrected solution on gird h can be written as 

Qh — Q'h + Vh (9-22) 


2h 


Q'Ah~ I Q'2h 

Ah 


(9-21) 


where 

k 

vh= T 

2ft 

This process is repeated until the solution converges. 


Q'2h- i Q'h 

2 h 


(9-23) 
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h (S,R 



S- Solve 
R- Restrict 
P- Prolongate 


Figure 9-1. Multi grid V-Cycle 

The sequence in which the transfers are done between grids can be done in a variety of ways. 
The usual strategy is to use the V-cycle or the W-cycle. In the V-cycle, the problem is solved on 
successively coarser grids, but driven by the finest grid residual, followed by the transfer of the 
correction and solution on successively finer grids all the way to the finest grid. Fig. 9- la shows a 
4-level V-cycle. In the W-cycle, intermediate V-cycles are performed on the coarser grids. Fig. 9- 
2b shows a 4-level W-cyc le. Computationally W-cycles are about 50% more expensive than 
V-cycle, however, W-cycles are considered more robust. The number of iterations of the solver 
required on each grid, N h jV 2 /,..., A/V - etc., to get the optimum performance, depends on the 
problem. In general the choice is to do more iterations on the coarser grids than on the finer ones. 


Restriction Operator 

The summation of the residuals of the fine grid cells, that make up a coarse grid cell is used 
as the restricted residual at that cell. Therefore, since there are 8 fine cells per course cell (Refer 
to Eq. 9-16), 

t [A -£*(<&)]= E [A-A(Oi)] (9-24) 

“ clfel 

The conserved variables are restricted by a volume weighted average. Thus 
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Figure 9-2. Multigrid W-Cycle 



E Q'h/Jh 

celt= 1 

E i /Jh 

ce//= 1 


(9-25) 


Prolongation Operator 

The corrections are prolongated to the fine grid cells using a tri-linear interpolation from the 
neighboring coarse grid cells. The interpolation function used to prolongate coarse grid cell-center 
values at A, B, C, D, E, F, G, and H to a fine grid cell-center ’a’, as shown in Fig. 9-3, can be 
written as 

h 27 9 3 i 

l h 8 Qu = g4 + ^Qb+SQd+SQe) + ^(SQc+SQ f + 6Q u ) + ±6Q Cl (9-26) 

where (Refer to Eq. 9-22), 


*<hh = Q'n- l Oi 

2h 


(9-27) 
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(^J) Coarse Grid Cell-center 


O Pi ne Grid Cell-center 


Figure 9-3. Coarse to Fine grid transfer 

Smoothing of Residuals and Corrections before transfer 

After restricting to a coarse grid, the residuals are smoothed by the implicit operator given 
below 


(l -e£<5^)(l -e^6 vv )( 1 -e<^)U = U (9-28) 

where %, S vr) and 6^ aje central- second-difference operators in the ij and k directions 
respectively. Thus, for example, 


- 2V i>j>k + 0i +u jt (9-29) 

U is the quantity that needs to be smoothed and 0 is the result of the smoothing. For example 
after restricting from grid h to 2 h the residual 

U=i [fh~ C h ((&)] (9-30) 

2 h 

(Refer lo Eq. 9-16) is smoothed. Similarly before prolongation from grid 2 h to h the corrections 

U = Q' 2h - i Q’ h (9-31) 

2/i 
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(Refer to Eq. 9-22) are smoothed 1361 . £{,e v and are free parameters that are used to control the 
amount of smoothing. Typically 


0 < £^£ v ,£(; < 1 


(9-32) 


Full Multigrid 

In order to obtain a good initial approximation for the finest grid, the problem, Eq. 9-12, 
is first iterated on a coarsest grid and the solution in prolongated to the next finer grid, several 
multigrid cycles are performed between these two levels and the solution is prolongated to the 
next finer grid. This process is repeated until the finest grid is reached where the solution is used 
as the initial guess and the multigrid process is continued. These preliminary multigrid cycles are 
inexpensive and results in a lower initial residual for the finest grid. This helps to drive the final 
residual to machine-zero in a fewer number of interations. 
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10. Flow Conditions 

Fluid Constants 


The fluid in this study is assumed to be air and the following constants are used. 

Prandtl Number, Pr = 1.0 

Ratio of specific he&ts, 7 = 1.4 

Exponent of Temperature to evaluate viscosity, ui = 1.0 


Incoming Flow Profiles 


The reference length used in the non-dimensionalization of the governing equations is the 
radius of the vortex-core at the incoming-face, x = 0. The Reynolds number is based on this 
reference length and a reference velocity = 1. The initial velocity profiles are the same one that has 
been used by several previous investigators including the present authors^ in their axisymmetric 
studies. They can be written in the cartesian coordinate system as follows 


U = Moo 

v = —MooVz/d 
w = MooVy/d 

where V is the velocity component in the 2 / 2 -plane, which is given as 

V=Sd(2-d 2 ), 0 < </ < 1 

V =S/d , d > 1 


( 10 - 1 ) 


( 10 - 2 ) 


where 

d= (y 2 + * 2 ) 1/2 (10-3) 

is the radial distance from the center of the vortex and S, the swirl parameter, is the circumferential 
velocity at the edge of the vortex-core. S is a measure of the strength of the incoming vortex. 
Like Moo and Re, it is a free parameter of the problem. 

In order to determine the profiles of the thermodynamic variables p and e a suitable total 
enthalpy was calculated based on the maximum total velocity and a temperature = 1. It can be 
shown that the maximum velocity occurs at d = (2/3 ) ly ^ 2 and the maximum velocity can be 
written as 


V*', = Ml(l + sh(2/3f) 

Therefore the total enthalpy can be written as 

H = l/o + Vlj 2 


(10-4) 


(10-5) 
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( 10 - 6 ) 



For Moo = 0.1 and S = 1, the profiles of some of the quantities, expressed above, are shown in 
the Figs 10-1 to 10-6. The quantities shown are plotted along y at (x,z : 0,0). The incoming-flow 
conditions are used as initial conditions at all the x locations. 



Figure 10-1. Density Profile Figure 10-2. Energy Profile 
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11. Results 


The cases that are presented here were run either on the Numerical Aerodynamic Simulator 
(NAS) Cray-2, called Navier, located at NASA Ames Research Center or on the Cray-2 called 
Voyager located at NASA Langley Research Center. The CPU time shown are normalized with 
respect to the Voyager CPU clock. 

All the solutions shown have a reference Mach no. =0.1. The computational domain, 
based on the vortex-core radius, is 10x8x8 in the x,y and z directions respectively and there 
are 64x32x32 cells in this domain with grid clustering around the vortex-axis. The boundaries 
around the vortex-axis (side-faces) are four vortex-core radii away from the x-axis. Fully upwind 
(k = -1) differencing with no limiting was used to discretize the inviscid fluxes. A three level 
V-cycle multigrid was used to accelerate the convergence. A global-minimum time step with a 
CFL = 5 was used to advance the solution in time. About 2500 multigrid cycles were required 
to drive the steady state residual to machine-zero. 


Case 1: Re = 100, S = 1 and = 0.1 


Fig 1 1-1 shows the residual history of this case run with and without multigrid. It is clear that, 
by using multigrid, an additional six order of magnitude reduction in steady state residual was 
obtained in the same amount of time compared to the single grid calculation. Fig 11-2 shows the 
breakdown of the vortex. A tube of particles seeded around the vortex axis, near the incoming- 
face, bulges around the breakdown region. A “tape” of particles seeded above the axis can be 
seen to follow the swirling flow, stretching and rolling-up as it travels downstream. It is clearly 
evident here that by using particle tubes and tapes, to visualize the flow, a lot more information 
about the flowfield can be obtained compared to tracking individual panicles. This solution 
agrees, qualitatively, with calculations done in the past^ 22 ^ using an incompressible axisymmetric 
formulation. Figs. 11-3 to 11-8 shows the profiles of some of the flow variables at the incoming- 
face along the y-direction at (x, z : 0, 0). These profiles are essentially axisymmetric. A slight 
asymmetry near the boundary can be noticed in the v-profile. This may be due to the boundary 
conditions on the side-faces. On the side faces, due to the swirl in the flowfield, one half of 
the face is treated as an inflow boundary while the other half is treated as an outflow boundary. 
Hence the entropy in the ghost cells on the inflow part of the face remains fixed at the initial level 
while the ghost ceils in the outflow part of the face get updated with entropy from the interior 
of the flow. If the side faces are sufficiently far away from the vortex-axis the effect this should 
be negligible. Fig 11-9 is a plot of the axial velocity along the x-axis. The region of negative 
velocities indicates the extend of the bubble. 
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CPU Hours (Cray-2) 

Figure ll-l. Residual History (Re = 100, 5 = 1, = 0.1) 



| 


Figure 11-2. Streamsurfaces (Re — 100, 5 = 1, Moo = 0.1) 
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Case 2: Re = 375, S = 1.464 and Moo = 0.1 

Fig. 11-10 shows the steady state residual history for this case. Fig. 11-11 shows the 
streamsurfaces. In this case the particle tube around the breakdown region has a mushroom like 
shape. The recirculation bubble has become pinched at the axis and is nearly shaped like a toroid. 
This pinching of the bubble has been observed in the axisymmetric calculations done in the past by 
the present authors^. A particle tape seeded above the axis shows a more pronounced stretching 
and roll-up as the particles move downstream. Fig. 11-12 shows the axial-velocity along the 
z-axis and, due to the pinching of the bubble, the flow along the axis does not have any negative 
velocities. However the axial-velocity along a line in the y-direction, at (x,z : 1.17,0), shown in 
Fig. 11-13, has negative velocities indicating a recirculation zone. 



CPU Hours (Cray-2) 

Figure 11-10. Residual History (Re = 375, S = 1.464, Moo = 0.1) 
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Case 3: Re = 400, 5 = 1 and M oo = 0.1 

Fig. 11-14 shows the steady state residual history for this case. Fig. 11-15 shows the 
streamsurface and particle traces. In this case a set of particles seeded near the axis can be seen 
to go around the breakdown bubble. A particle seeded inside the bubble stays trapped inside 
the recirculation region. Again the observations are consistent with the ones seen in previous 
axisymmetric simulations. A particle tape seeded above the axis shows stretching and roll-up as 
the particles move downstream. 



CPU Hours (Cray-2) 

Figure 11-14. Residual History (Re = 400, 5 = 1, Moo = 0.1) 



Figure 11-15. Particle Traces/Streamsurface (Re = 400, 5 = 1, Moo — 0.1) 
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12. Concluding Remarks 

A three dimensional iinite volume full Navier-Stokes code has been developed to study the 
vortex breakdown phenomenon. The code is robust and its performance, both in terms of computer 
resources requirement and rate of convergence, is good. 

The bubble-type breakdown has been simulated for several sets of Reynolds numbers and 
Swirl velocity parameters. The results agree, qualitatively, with the ones obtained in the past 
using an axisymmetric incompressible formulation. We believe, that this code can be used as a 
good platform to conduct further extensive numerical experiments and to investigate other types 
of breakdown. 
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Appendix A 


Consider a cell i,j,k in the flowfield, shown in Fig. 3-3. Let Q t j k represent the conserved 
variables located at the cell-center. The semidiscrete finite-volume representation of the governing 
equations is (Eq. 3-4) 


-<*-**- iJ* 

+( F -n, j+ ik-(F-n, H>k 

+(G - G)i >Jik +i - (G - G) hhk _ i = 0 


(A-l) 


where the fluxes are evaluated at the cell interfaces defined by the grid points. Substituting Eqs. 2- 
24 and 2-25 we get 




f .ISAM 


+ 


+ 


fa v^iv^li 


f VOIZO 
' ivci) J 
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F ■ 




- '( F '' 


vn 

|Vf|l 

m) 

J J 

Vr, ' 

IVvIl 


1 J 


liven 

|V(|> 

1 J J 






= 0 


(A-2) 


j 


where 


f = (e - e)i + (f - fp + (g- g)k 


(A-3) 


It can be shown that the terms £a;/|Vf|, £y/|V£|, £*/|V£| are the components of the unit normal 
vector, pointing in the ^-direction, from the cell face that is normal to the ^-direction. |V£|/J is 
the directed area of the cell-face and 1 jJ is the volume bounded by the cell-interfaces. Hence 
Eq. A-2 can be written as 
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I $2' i + 

J * >i,„k 


+ 


^■“) |AS| ],„ +! , i + K fs ) |4S| ] i ,H. t 


(A-4) 


where 


h = n x i + n y j + n z k (A-5) 

is the unit normal vector pointing outward from the cell-face and |A5| is the cell-face area. If the 
conserved variables Qij,k is regarded as a cell-averaged value, i.e., 




y J J J Q(x,y,z,nAt) iJik dV 


(A-6) 


where V is the cell volume, Eq. A-4 can be written in the integral form as 


ltIJi QdV+ JJ pAdS=0 (A ' 7) 

v s 

This is the integral form of the conservation equations. Hence Eq. 3-4 is a consistent approximation 
to the integral form of the conservation equations. 
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Appendix B 


The matrices Trll+T^ 1 and Trll'T^ 1 are the same as the matrix T^T;: 1 except that 
they are evaluated using Roe-averaged quantities and the eigenvalue matrix contain only either the 
positive or negative eigenvalues. In order to evaluate the matrix T r A r T^ ! we need to know the 
tangential components of the velocity v and w that are parallel the unit vectors l and m. However 
the choice of 1 and fh is arbitrary. The arbitrariness in choosing the tangential vectors can be 
avoided in the following way. By definition we can write 

f_x Ty f z u 

l x ly l Z v (B-l) 

fh x fh y fh z w 

Since r,l and m are mutually perpendicular unit vectors the 3x3 matrix in Eq. B-l is an orthogonal 
matrix. Hence its inverse is the same as its transpose, i.e., 

r_x r y f_ z ' 

L ly T z = I (B-2) 

_m x fhy fh z 

Therefore we can write 

u f x l x fh x u 

V = fy ly fhy V (B‘3) 

w f z l z fh z w 

From Eq. B-3 we get 

vl x + wm x = u — uf x 

vly + Wfhy — V — ufy (B'4) 

vl z + wm z — w — uf z 

Eq. B-4 

Using Eqs. B-2 and B-4 and from the fact that 

f (B-5) 

the columns of matrix TrArT^ 1 can be written as follows. 
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Column 1 = 


oq 2 X T \ — u \ r 2 + A r s 
{oq 2 H - a 2 u 2 } A r i + u{oq 2 - H) A r2 
\(jq 2 u - a 2 ur x }X T i + {oq 2 f x - uu}X T 2 
\oq 2 v - a 2 uf y }X r i + \oq 2 f y - uu}A r 2 
{oq 2 w - a 2 uf z ) A r i + {crq 2 r z ~ uw}X r 2 , 


Column 2 = 


aA r i 

oHX r \ + cruX r 2 + A r s 
ouX T \ + or x X r 2 
ovX T \ + of y X r 2 
owX r \ + 0r z X r 2 


Column 3 


— ouX r \ + fx X r 2 

{a 2 uf x - ouH] A r i + {Hr x - ouu}X r 2 
{a 2 f x f x — cnxu}A r i + {r x u — our x }X r 2 + A r ,i 
{a 2 r x f y — crm>}A r i + {f x v — our y ) X r 2 
{ a 2 f x r z — euw}X T i + {r x w — our z }X T 2 


Column 4 = 


— crvA r i + r y A r 2 

{ a 2 uf y - ovH} A r i + { Hf y - ovu}X r 2 
{ a 2 f y f x - ovu)X r \ + {f y u - ovr x }X r 2 
{a 2 f y r y — ovv}X T \ + {f y v — crvf y ) A r 2 + A r 3 
{a 2 f y f z - auu>}A r i + { f y w - ovr z }X r 2 


(B-6) 


where 


Column 5 = 


-cru>A r i + r z X t 2 

{< a 2 uf z - owH) A r i + {Hf z - owu}X r 2 
{ a 2 r z f x — owu) A r i + { f z u — crwr x }X r 2 
{a 2 f z f y - owv) A r i + {r z v - owr y }X r 2 
{a 2 f z f z - <m>u>}A r i + {f z w - owf z }X r 2 + A r s 


A r i = (A| + A7-2A° r )/2a 2 
Xr2 = (A+ - X-)j2a 
A r 3 = A® 


(B-7) 


Note that by substituting for A+,A; and A° r from Eq. 5-47, we can recover the matrix K r given 
by Eq. 5-39. 

From Eq. B-6 the matrices Trll+T;: 1 and Trl^T,: 1 can be evaluated by replacing the flow 
variables by their Roe-averages and by setting the negative and positive eigenvalues in II r to zero 
respectively. The evaluation of 11+ and II7 can be automated by using the following equations 
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(B-8) 


n+ = 

nr = 


n r + |n r | 
2 

n r - jgr| 
2 
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Appendix C 


Carrying out the linearization procedure outlined in section 6, we can write the viscous flux 
Jacobian matrix as 


where 



1 

o 

0 

0 0 O' 

v r J Vr ' M % 

r J Re ^ 

*21r 

*22r 

*23r *24r *2fir 

^31r 

0 

®33r ^34r ^35r 

*41r 

0 

*43r ^ 44r tf 4Sr 


. *51r 

0 

^53r ^54r ^55r. 


*31r 

’J'a.tr 

^34r 

^35r 


U V 

<Alrr ^2rr — 

P P 

4 > lrr 


P 

< hrr 

P 

4>?>tt 


<f> 3rr 


W 

P 


'J^lr — ^4 tt fibrr ^6rr 

P P 


W 

P 


^43r 
^ 44r 
^45r 


^4rr 

P 

<H rr 

P 

<Hrr 

P 


^51r = ~hrT~ 

P 


^53r = 


^7rr 

P 


$54r = 


P 


^55r = 


<farr 

P 


<f>Zrr 


V 

P 


farr 


W 

P 


(C-l) 


(C-2) 


(C-3) 


(C-4) 
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*21 r = *31r^ + * 41r U + *51rW ~ 


<^0rr7 ( e_2 9 2 ) 


*22r 


<ftOrr7 1 

P r p 


*23r = *33r“ + *43r*> + *53r^ “ 
*24r = *34r« + *44rl> + *54rW - 


*25r = *35rU + *45rU + *55rU> ~ 


^Orr7 u 
Pr P 
farrl V 
Pr P 

hrpyw 

Pr P 


and r stands for £, r) or (. Using Eqs. 2-25, 6-7, 6-8, 6-9 and 6-10 in Eqs. 
we can write 


*31r = J 


|Vr 


|Vr| 


*34r = 


_ r|Vr|l/l 

5r “ 3 J p \ 3 
| Vr 


*35 
*41r = J 


*43r = J 
*44r = J 
*45r = J 

*51r = J' 
*53r = J 
^ 54r = J 
*55 r = J 


J 

|Vr 

J 

|Vr 

J 

|Vr 


J 

| Vr 


|V» 


J 

|Vr 


J 

|Vr 


If 1 . _ _ J 

-< — u - -r x (ur x + vr y + wr z ) j 


1 / 1 - - 1 

~pW x> 


} 


- ^r y (ur x + vr y + u;r z ) j 

HM 

^ j-tu - if y (uf* + vr y + urr 2 ) j 


Mr 

H i+ i ,A } 


-r 2 r v 


(C-5) 


C-2, C-3, C-4 and C-5 


(C-6) 


(C-7) 


(C-8) 
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where 


*21r 

^22r 

#23r 

*24r 

^25r 


J ~-^-2q 2 - ^(uf^ + vr y +wr z ) 2 - ^(e-2 q 1 

t Vr l i r_7_\ 

J p \ Pr / 

_ jVrl If 1 . _ _ . 7 ^] 

J — + gMur* + wr„ + wr 2 ) - — j 

_ IVrl if 1 _ . _ _ _ N 7^1 

J — -|v + -r v (ur* + vr y + wr z ) - — | 

T Vrl If 1 _ . _ _ . 7U>1 

J -y-- ju; + -r 2 (ur x + ur y + u>r 2 ) - — J 


(C-9) 


9 


2 


u 2 + 1> 2 + u> 2 
2 


(C-10) 


Substituting Eqs. C-6, C-7, C-8 and C-9 in Eq. C-l we get the final form given in Eq. 6-25. 
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